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ON THE ANALYSIS AND CONSTRUCTION OF PERFECTLY MATCHED LAYERS 
FOR THE LINEARIZED EULER EQUATIONS 


J. S. HESTHAVEN * 

Abstract. We present a detailed analysis of a recently proposed perfectly matched layer (PML) method 
for the absorption of acoustic waves. The split set of equations is shown to be only weakly well- posed, 
and ill-posed under small, low order, perturbations. This analysis provides an explanation for the stability 
problems associated with the spilt field formulation and illustrates why applying a filter has a stabilizing 
effect. 

Utilizing recent results obtained within the context of electromagnetics, we develop strongly well- posed 
absorbing layers for the linearized Euler equations. The schemes are shown to be perfectly absorbing in- 
dependent of frequency and the angle of incidence of the wave in the case of a quiescent mean flow. In 
the general case of a convecting mean flow, a number of techniques are combined to obtain an absorbing 
layer exhibiting PML-like behavior. The efficiency of the proposed absorbing layer is illustrated through the 
computation of benchmark problems in aero- acoustics. 

Key words, acoustics, absorbing boundary conditions, hyperbolic systems 

Subject classification. Applied and Numerical Mathematics, Fluid Mechanics 

1. Introduction. When addressing wave- dominated problems, as found in aero- acoustics or electro- 
magnetics, one is often facing the problem of how to accurately obtain infinite domain solutions using a finite 
computational domain. The truncation of the computational domain must be done in a way that avoids, at 
least approximately, the excitation of reflected waves which might otherwise contaminate the computational 
domain and falsify the solution. 

The issue of how to properly devise such boundary conditions at an artificial computational boundary 
has received much attention in past. The use of characteristic boundary conditions [1] is appealing due to 
its simplicity, but is only accurate for close to perpendicular incident of waves. More elaborate schemes 
involve radiation boundary conditions based on localization of the Dirichlet-to-Neumann map [2, 3] or an 
asymptotic expansion of the far-field solution [4]. A fairly recent review of such methods can be found in [5]. 
Alternatives to such schemes involve the introduction of buffer or sponge layers in which the waves are either 
damped [6], accelerated to supersonic conditions [7], decelerated [8] or attenuated by combinations thereof 
[9]. The construction of such schemes are in most cases based on physical arguments with little theoretical 
evidence of their, often quite remarkable, performance. 

In the context of electromagnetics, Berenger [10] recently proposed a novel way by which to derive the 
sought after absorbing boundary conditions. By splitting Maxwell’s equations in an unphysical manner, ad- 
ditional degrees of freedom are introduced thus allowing for the construction of absorbing layers. These layers 
have the remarkable property that they maintain their absorbing properties irrespective of the frequency 
and angle of propagation of the incident wave, i.e. this approach appears to provide an optimal absorbing 
boundary condition. These layers, called Perfectly Matched Layers (PML), are now under intensive research. 
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research was partially supported by the National Aeronautics and Space Administration under NASA Contract Nos. NAS1- 
97046 and NAS 1-19480 while the author was in residence at the Institute for Computer Applications in Science and Engineer- 
ing(ICASE), M/S 403, NASA Langley Research Center, Hampton, VA, 23681-0001 
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While the PML scheme has been applied successfully during the last years, it was recently proven [11] that 
the splitting of Maxwell’s equations makes the resulting set of equations only weakly well-posed and ill-posed 
under arbitrary low order perturbations, i.e. the numerical schemes resulting from these equations can be 
expected to be unconditionally unstable, an example of such being provided in [11]. This realization has 
focused the attention towards alternative well- posed formulations of the electromagnetic PML methods and 
several such schemes have been proposed in recent years, see. e.g. [12, 13, 14, 15, 16, 17]. Hence, although 
the original approach for the construction of PML schemes have proven erroneous, the general approach have 
proven extremely fruitful and allowed for the computation of problems in electromagnetics of unsurpassed 
accuracy. 

Inspired by the success of the PML methods for Maxwell’s equations, Hu [18] recently proposed a PML 
method for the equations of acoustics by taking an approach very similar to the one originally developed for 
Maxwell’s equations, i.e. by splitting the equations of acoustics in an unphysical manner. However, contrary 
to most work within the community of electromagnetics, Hu [18] reported the need for using a low pass filter 
inside the absorbing layers to maintain stability of the scheme. A similar observation was made in [19] where 
no filter is applied and the numerical solutions are found to exhibit exponential growth. This points to an 
inherent instability of the scheme for which a partial explanation, in terms of loss of strong well-posedness 
in certain special cases, is providedin [19]. 

In the present work we provide a complete analysis of the split PML scheme of [18], confirming the 
speculations put forward in [19] in a more general context. Indeed, the scheme of [18] is found to be only 
weakly well-posed in the two-dimensional case and ill-posed under low order perturbations. We proceed by 
presenting a well-posed PML scheme for the quiescent equations of acoustics, and a well- posed absorbing 
layer exhibiting PMI^like behavior for the general convective case. 

The remainder of this paper is organized as follows. In Sec. 2 we introduce the equations of acoustics 
as obtained from the linearized Euler equations. Section 3 contains the first part of the paper in which we 
present an analysis of the PML method recently proposed in [18] and provide an explanation for the problems 
of stability reported in [18, 19]. This leads to Sec. 4, where we present an alternative to the unstable PML 
scheme. For the case of a quiescent mean flow we construct a well behaved PML method and illustrate its 
performance through numerical experiments. For the general case of a convective mean flow, we propose 
to apply a combination of techniques to arrive at absorbing layers with PML- like behavior and support the 
reasoning by numerical studies. Section 5 contains a few concluding remarks. 

2. The Equations of Acoustics. We shall consider the two-dimensional compressible Euler equations, 
linearized around a uniform parallel mean flow of the form (poi ^o, 0,po)- Choosing v$ = 0 does not introduce 
any restrictions on the analysis as the general situation may always be rotated to arrive at this particular 
case. Within this scenario, the equations take the form 


(i) 


dq .dq _ dq 

m +A di +B di ~ 0 ’ 


where the state vector, q, and the constant matrices, A and B, are given as 
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These equations are recovered from the Euler equations by linearizing around the uniform mean state, 
(po,uo 5 0,po), and introducing the normalization 
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where L represents a characteristic length while Co = yJPo/po refers to the sound speed of the mean flow. 
In this context, Eqs.(l)-(2), describes the propagation and interaction of waves in a parallel uniform flow 
with a Mach number, M = uq/cq. 

A deeper understanding of the underlying properties, physical as well as mathematical, of Eqs.(l)-(2) 
can be gained by introducing the similarity transform 
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where S~ 1 q = R = [p + u, p — p, v,p — u] T represents the characteristic variables. We recognize the con- 
vective entropy ( R<i ) and vorticity waves (R 3 ), respectively, and the co- ( R \ ) and counter- propagating ( R 4 ) 
sounds waves through which the complete physical scenario can be understood. 

However, the similarity transformation also shows that A and B can be transformed such as to become 
symmetric simultaneously by multiplication with a positive definite diagonal matrix. Within the context of 
the present work this has the consequence that Eqs.(l)-(2) forms a strongly well-posed hyperbolic system 
[20] implying that the well-posedness of Eqs.(l)-(2) is unaffected by the addition of low order terms. As we 
shall see shortly, lack of strong well-posedness can have serious consequences and make the construction of 
convergent numerical schemes impossible due to inherent instabilities of the system of equations. 


3. An Analysis of the Split-Field PML Method. Following the line of thought initiated in [10] for 
the development of perfectly matched layers (PML) for electromagnetics, Hu [18] recently proposed a split- 
field PML scheme for the two-dimensional linearized Euler equations as given in Eqs.(l)-(2). Different from 
the approach of [10], in which only some of the field components are split, in [18] all the field components of 
q are split to arrive at a set of equations to be solved in the layer of the form 


( 3 ) 
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+ B® + C 3 q s = 0 
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where q s = [pi, p 2 , U\,U 2 , t>i, t» 2 ,Pi,P 2 ] T ? such that p = p\ + p 2 and likewise for the velocity components and 
the pressure. The 8x8 matrices in Eq.(3) become [18] 
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while C s = diag(or x1 Oy,a x ,<7 x ,<jy,a x ,cr x ,<jy) represents the diagonal matrix responsible for the dissipation 
of the waves. 

At first the use of split variables may seem perfectly legal since for a x — a y = 0 the original equations 
are recovered by adding the equations for the split fields. This reasoning, however, is faulty as we shall show 
in the following. 

Let us first address the issue of well-posedness of the split system of equations, Eqs.(3)-(4) and recall that 
the question of well-posedness of the system is unaffected by the low order term, C s q 3 , which we therefore 
neglect in the subsequent analysis. We first consider the diagonalizing similarity transform of A 5 given as 
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to obtain S -1 A S S = diag (0, 0, 0, 0, M — 1,M, M, M + 1), i.e., 4 zero eigenvalues have been introduced as 
a consequence of the splitting. Now, if S and S -1 cannot transform B 5 into a matrix that can be made 
symmetric by multiplication with a positive definite diagonal matrix, then A s and B 5 cannot be symmetrized 
simultaneously [21]. Indeed, we obtain 


0 0 0 0 0 1 0 0 

0 0 0 0 -1 0 0 -1 

0 0 0 0 0 0 0 
0 0 0 0 0 1 0 0 

00000 \ 00 

0 0 0 0 1 0 0 1 

00000 0 00 

00000 \ 00 


which, due to the zeros in the left half, clearly can not be made symmetric. This observation, however, is not 
conclusive in terms of lack of strong well-posedness of Eqs.(3)-(4), but it remains a bad sign, since the split 
set of equations has lost an important symmetry property as compared to the original un- split linearized 
Euler equations, Eqs.(l)-(2). 
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To continue the analysis we shall focus the attention on the Cauchy problem, i.e., neglect the effect of 
the boundary conditions in Eqs.(3)-(4). We introduce the spatial Fourier transform of q s as 

/ oo poo 

/ dk x dky , 

-OO J — OO 

where q s = [pi, p 2 , &i, U 2 ? vu V 2 ,Pi,p 2 } T represents the Fourier coefficients of the split field components to 
arrive at the initial value problem 


(5) 
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where the symbol, P (ft x , k y ), is given as 
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Integration of Eqs.(5)-(6) is achieved by first realizing that the evolution of the individual split components 
depends only on the un-split variables. Hence, by obtaining the solution of the Cauchy problem of Eqs.(l)-(2), 
given on the form 
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we may obtain the solution for the split field variables by introducing the solutions of the un- split variables 
into Eqs.(5)-(6). 

Considering the initial conditions 


9(0) = [po,w 0 ,t)o,Po] T , 


we obtain the solution to Eq.(7) as 
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with the three vectors a = [a p , a u , a v , a p ] T , b~ [b p , b UJ b v , b p } T and c= [c p , c^, Ct,, Cp] T , having the entries 
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Here we have 


v = 



(i = u 0 k x + v 0 k y , 


and we have utilized that po = po for consistency. In Eq.(8) we immediately recognize the three types of 
waves, inherent in the linearized Euler equations, giving rise to three different wave speeds. Moreover, we 
note that a and b as well as c are bounded for all values of k x and k y confirming that Eqs.(l)-(2) constitute 
a strongly well-posed problem for which the solution can be bounded up to exponential growth in time by 
the norm of the initial data. 

Let us now continue in order to arrive at the solution of the Cauchy problem for q 2 = [p 2 , u 2 ,V 2 ,P 2 ] T > 
bearing in mind that we could equally well have chosen the other set of equations. The dynamics of these 4 
variables are described by the system 
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Integrating the solution, Eq.(8), we recover the solution to Eq.(9) on the form 
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For the spht set of equations to be strongly well-posed we must ensure that the solution, Eq.(10), is bounded 
by the norm of the initial data for any choice of k x and k y , or, in other words, a£, b 2 and must remain 
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bounded for any combination of k x and k y . It is easily verified that and indeed remain bounded for all 
values of k x and k y provided the mean flow is purely subsonic, i.e. \M\ < 1. However, absorbing boundary 
conditions are not necessary in the case of supersonic flow as reflections from the open boundary are unable 
to enter the computational domain and therefore will not cause any problems. 

The situation for is very different. In the limit where \Mk x \ — > 0 and \k y \ \Mk x \ we can only 
bound two of the terms in as 


( 11 ) 




< \Cvky\t , 


i.e, we recover terms that grow in time with a coefficient, k yi which is unbounded. Hence, [|p2 1| and ||p2|| 
cannot be bounded by the norm of the initial conditions, but rather depends also on the norm of the 
derivatives of the initial conditions, iio and t)o. Consequently, the split set of equations, Eq.(3), is only 
weakly well-posed with the solution depending not only on the initial conditions, but also on the smoothness 
of these data. 

It is noteworthy that, as is the case for the split field perfectly matched layer methods of electromagnetics 
[10, 11], in the case where |Ar y | = 0 strong well-posedness is recovered. Hence, the one-dimensional version 
of the split field method for the perfectly matched layers of acoustics is valid and well suited for numerical 
solution. 

The loss of derivatives is as such not a severe problem. However, contrary to strongly well-posed 
hyperbolic problems, it is well known that weakly well- posed systems may become ill-posed under low order 
perturbations [20] , thus rendering the systems of equations inherently unstable and proper numerical solution 
impossible. 

To see this, we introduce a low order perturbation of the form 
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q s 


i.e. the perturbation corresponds to a small perturbation in the split field velocity component, V\ and t? 2 , 
however maintaining that v = v\ + • 

We consider the perturbed Cauchy problem 


dt 


(P(/c x , k y ) + E) q s =Pq s 


A necessary condition for the perturbed problem to remain well-posed is that the real parts of the eigenvalues 
of P remain in the left half plane for any choice of k x and k y and, preferably, also for any choice of £ and 
M. The 4 eigenvalues of P are given as 


Ai = A 2 — A3 = 0 , A 4 — -ik x M 
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while the remaining 4 eigenvalues appear as the roots of a complex polynomial 


/4(A) = O4A 4 + O3A 3 ■+■ Q2 A^ -f- OqA + C*0 = 0 ) 
in which the coefficients are given as 

0:4 — 1 , Q3 = e + i3MA: x , c*2 = * 2 + * 2 - 3M 2 * 2 + ie(k y + 2k x M) , 


<*1 = e(k 2 x + _ 3Mh x k y - M 2 * 2 ) + iM* x (* 2 - M 2 * 2 + * 2 ) , 


a 0 = i2eMk x ky(l — M) . 

Rather than solving the complex polynomial to obtain the expressions for the eigenvalues, we shall recall 
the Routh-Hurwitz criteria expressing that all the roots, A, of an n’th-order polynomial, P n ( A), lies in the 
left half plane if and only if the roots of the (n — l)’th-order polynomial 

n 

Pn-i(A) = [a n a* n _ 1 +a n - 1 a* n - a„<A] P„(A) + A ^(-1)"“^ , 

1=0 

lies in the left half plane and the real part of a„_i/a n is positive. Hence, by successively applying this result 
we arrive at criteria under which the perturbed initial value problem remains well-posed. 

The first condition immediately yields the requirement that Re(a3/a4) = e > 0 which certainly is a bad 
sign since we can not in general control the sign of the perturbation. However, one additional application of 
the Routh-Hurwitz criteria results in a condition for well-posedness as 

\Mk x \ >1*^1 . 

This condition is very similar to the limit for boundedness of C2, Eq.(ll), and confirms that the weakly 
well-posed system, arrived at by splitting the linearized Euler equations to facilitate the development of 
the perfectly matched layers, becomes ill- posed under low order perturbations. We should note that there 
is nothing particular about the low order perturbation, E. Indeed, ill-posedness can be shown, using the 
technique outlined above, for perturbations of the velocity components as well as of the density and pressure 
components. 

Due to finite precision, low order perturbations will always exist in actual implementations of the split 
field equations, i.e., problems with maintaining stability of these schemes should be expected. Indeed, this 
is exactly what was reported in [18, 19] where it was found that applying a filter in the PML layers was 
necessary to maintain stability. An indication of why the filter has a stabilizing effect for this problem is 
provided by the condition for boundedness, Eq.(ll). Indeed, if the filter is sufficiently strong as to ensure 
that \Mk x \ > Ifcyl for ail values of |* x | and |* y |, e.g. by enforcing a strong filter along t/, the system remains 
well-posed and, as a consequence, the scheme might recover stability or at least postpone the effects of the 
instability. 

4. The Construction of Well-posed PML Methods. The inherent instability appearing as a result 
of the splitting of the linearized Euler equations illustrates the care that has to be exercised in devising 
absorbing layers for such types of equations. 


8 



The weakly well-posedness and associated ill-posedness under small perturbations of the split-field PML 
equations was recently shown [11] also for the original PML method as proposed in [10] and several successful 
attempts have been made to formulate strongly well- posed PML methods for the equations of electromag- 
netics [12, 13, 14, 16, 17]. Hence, rather than attempting an ab initio development of perfectly absorbing 
layers for the linearized Euler equations, we shall utilize this recent development to arrive at the sought after 
well behaved methods. 

A strongly well-posed PML method for Maxwell’s equations is proposed in [17] and tested numerically 
in [16] and we shall base the remaining part of this paper on this particular formulation. We should em- 
phasize though that alternative well-posed formulations might equally well be employed as the basis of the 
development of the PML methods for the equations of acoustics. 

In what remains we restrict the attention to the case of a purely subsonic free stream, i.e. \M\ <1 
in Eq.(l). This poses no restrictions on the analysis and the applicability of the proposed schemes as the 
situation for \M\ > 1 is trivial since all information is leaving the computational domain and, as information is 
prohibited from propagating upstream, reflections cannot enter the domain and interfere with the computed 
results. 


4.1. The Quiescent Case. Let us first consider the simple case of a quiescent free-stream, i.e. M = 0, 
for which the equations are given as 

du dv 
dx dy ’ 
dp 

dx 1 
dp 

dy 1 
du dv 
dx dy 

appearing directly from Eqs.(l)-(2). 

We propose to consider an absorbing layer for the quiescent linearized Euler equations of the form 
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Here, e — e(x) and p — p(y) signifies the non-dimensional damping parameters. We immediately note 
that since the Euler equations are altered by the introduction of low-order terms only, the system of partial 
differential equations is well-posed by construction while the additional freedom required for obtaining the 
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sought after properties of the matched layers are introduced through 4 additional equations describing the 
development of the artificial fields, P x and Q x , along x and, likewise, P y and Q y) along y. In general, we 
assume that the absorbing region is outside a square bounded by |x| = a and \y\ = b while the specification 
of the parameters, e and p remains open at this point in time. 

To come to an understanding of the properties of this absorbing layer, we follow the analysis introduced 
in [11] and study the behavior of a plane wave hitting the layer interface, which we assume is positioned 
at x = 0. As the system we are dealing with is purely linear, it poses no restrictions to consider only the 
behavior of plane waves, since any type of excitation can be decomposed into a superposition of such plane 
waves of the form 


(14) 
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where c? + /3 2 = 1 represents the arbitrary angle of incidence and u> signifies the normalized frequency of 
the incoming wave. 

We shall seek solutions inside a layer in the x-direction, i.e p(y) = 0 in Eq.(13), of the form 


(15) 
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Qx 


Qx(x) 

Px 


_ Px(x) _ 


and shall obtain the full solutions for the fields inside the layer. 
Introducing Eq.(15) into Eq.(13) yields the equations 


(16) 


1 £ 

P = P , V = PP , Q X = -T« , Px = — w , 

y IU) 


thus expressing 4 of the 6 variables in terms of u and p . The latter are governed by the coupled equations 

d ~ i> 2 - d ■ , 2-. 

—p = — — u , — \ipu\ = —luipa p , 
ax iu) ax 

where we have introduced ip(x) — iuj 4- e(x). Combining these two equations yields a second order variable 
coefficient ODE for p on the form 


A 

dx 


a± P ) 

dx F J 


ifjctp , 


which has an analytical solution given by 


(17) p(x) = Ae° So *•»> d " + Be~ a SS dr > , 

through which we immediately obtain the solution to u as 
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( 18 ) fi(ar) = -~^~p{x) = -a^ ( Ae a £ dT > - Be ~ a £ *A . 

\p z ax W \ ) 

The remaining fields are then given from this using Eq.( 16 ). 

The specification A and B naturally depends on the boundary conditions we choose to impose and, 

indeed, there are several ways of doing so. We shall assume that the layer has a finite width, d, and shall 

hence need to impose boundary conditions at x — 0 and x — d. For the solution of hyperbolic systems it is 

most natural to impose characteristic boundary conditions by specifying the incoming characteristics. This 

amounts to requiring that R\ = p+u remains continuous across the interface x = 0 while R4 = p — u = 0 at 

x — d, i.e. no information is entering the layer. Imposing these boundary conditions, using Eq.( 14 ), implies 

p(0) + u(0) = 1 + q , p(d) — u(d) — 0 , 

from which we arrive at (a / 0) 


( 19 ) 

where 


A = - — -e~ 2a! B 


B = 


OL + 7 


1 1 a-7 1 -q -2 al ’ 

^ q +7 1 -fa 


7 = 


c(d) -+* \ 


, /== / dr} — iujd + / £{rf) dr] . 

do do 


Combining Eq.( 16 ) with Eqs.( 17 )-( 19 ) yields the complete solution inside the layer, 


( 20 ) p{x,y,t) = B 


l 2L 2 e 2t ^«(^-d) e “ 2a e{v)d*l 

0 + 7 


u{x,y, t) = —Bo. 
v{x,y,t) = (3B 
p(x,y,t) = B 


iuj 


^ __ 2L 2 e 2iw «( x - d ) e _2Q f x 

e(x) -+■ iuj |_ 0 + 7 

1 a ~ 7 ^2ia;Q(x— <q^~ 2t * /J* 
a + 7 

1 4 . Q ~~ / 7 c 2iwa(x-d) c -2Q /%(r?)dT) 

a + 7 


tw(t-ax-/3y) e - a f Q e(v) dr, 

e iu{t-<xx-Py) e ~OL f* e(T)) dr, 
e iu)(t-OLX-0y) e ~<x f Q £(r?)dV 
iu(t-ax-py) e ~<* f Q e(r,)drj 


P x (x,y,t) = -Ba 


e(rr) + iu; [ a + 7 


2 __ ^ 7 e 2 iwa(x-d) e - 2 « / x d £(^?) 


Qx(^ ? y,^) = 5 a- 


1 ® 7 e 2iufa ( 3:— ^) e _2a / x ^ 


e iv(t-ax-0y) e -<xf*e(v)<in 

iu>(t-cxx-py) g— a f Q e(v) dr. 


Since 


c-t-ia; 


(e(:r) + zu;) 1 2 [ a + 7 

< 1 all components in the layer are bounded by |p| which is bounded like < 0 with 

equality only for gazing waves, i.e. a — 0. Thus, all waves are damped independent of frequency and angle 
of incidence as should be required by a truly perfectly matched layer. 

We observe that the fields at the layer interface, x = 0, in general are discontinuous with a jump 
proportional 

1 4. Q ~"r c -2iujad -2al 
a +7 


1 _L_ a-7 X-q -2 al 

x ^ a +7 1 +a c 


which, however, is exponentially small. 

Naturally, an analysis equivalent to the above can be completed for a PML layer in the y-direction while 
a corner region, in which e > 0 as well as p > 0, can be analyzed using separation of variables, yielding 
results similar to those obtained in the layer above. 
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4.1.1. A Numerical Example. In order to confirm the theoretical analysis put forward in the previous 
section and study the efficiency of this new PML method, we have implemented the scheme on an equidistant 
grid using a 4’th order, centered, finite- difference scheme with 3’rd order closure for stability in space, while 
we use a 4’th order Runge-Kutta scheme for advancing the system, Eq.(12), in time. The time step, At, is 
chosen to be well below the stability limit. 

We note that contrary to the scheme proposed in [18], there is no need for applying a filter to maintain 
stability and, to emphasize this point, we have not used any filters in the present work. 

The initial conditions are taken from a benchmark problem of computational aeroacoustics found in [22] , 
namely 


(21) 


-(In 2) 


p(x,y) = e 
u(x,y) = 0.05 (y-yb)e 
v(x, y) = — 0.05(x — Xb)e 
p{x,y) = e 


(i-i») 2 -Kv- t q) 2 -(In 2) 

+ O.le 

-(in 2) l*-* 


(x-x b ) +(y-y b ) 

2 


~(ln2) (x ~ g fr )2 t (V ~ Vb)2 


(ln2) (x ~ T ^ a V V ~ Vaj2 


where ( x a , y a ) signifies the center of the initial sound pulse of width 6 a , while (x&, yb) refers to the center of 
the initial vorticity and entropy pulse of width The analytic solution to this problem may be obtained by 
exploiting the axis-symmetry of the initial conditions together with the use of Fourier transformation. The 
exact expressions for the solution are given in [22]. 

The profiles, e(x) and /i(y), required in the specification of the scheme, Eq.(13), are taken to be of the 
form 


(22) £(*) = C x — -) , p(y) = C y . 

\ XPML ) \ VPUL ) 

Here we have assumed that the computational domain in which Eq.(12) is solved is bounded by |x| < a 
and |y| < b while xpml and ypML refers to the width of the absorbing layers along x and y, respectively. 
The constants, C7 X , C y and n, controls the strength of the layer and we have chosen these parameters as 
C x — C y = 2 and n = 3. The auxiliary equations of Eq.(13) are advanced in time using the same scheme 
and time-step as for Eq.(12). 

We consider the problem in the computational domain (x, y) G [— 50, 50] 2 with the absorbing layers 
outside and position the acoustics pulse at (x a ,y a ) = (—25,0) with a width of S a = 3 while the non- 
propagating vorticity /entropy pulse is positioned at (x&, yb) = (25, 0) with a width of 6b = 4. The absorbing 
layers are terminated using characteristic boundary conditions as discussed during the analysis of the scheme. 

In Fig. 1 we show the pressure field at various times as computed using Ax = A y = 1 and At = 1 and 
xpml = ypML — 10, i.e. 10 computational cells in the absorbing layer. 

As expected from the analysis, the sound wave propagates undisturbed out of the computational domain 
with no visible reflections. The high frequency noise visible on the contours is a consequence of the accuracy 
of the scheme and the lack of filtering, rather than a result of reflections as can also be observed on Fig. 2, 
where we show the u - velocity field propagating undisturbed out of the computational domain. 

To verify the dependency of the efficiency of the absorbing layer on the width of the layer, we have 
computed the maximum pressure error along the line x = —48 as a function of time. In Fig. 3 we show 


12 



a) 


b) 




c) d) 




Fig. 1. Pressure contours for a sound pulse , propagating in a quiescent medium. The contour levels are ±0.1, ±0.05, ±0.01 
and ±0.005 with the computed result given at t — 20 (a), t = 40 (b), t = 60 (c) and t = 80 (d). 

the development of the pressure error for various layer widths as compared with using only characteristic 
boundary conditions to terminate the computational domain. 

Indeed, as expected, we see that even for a layer of only 6 cells the PML scheme out-perform the 
characteristic BC in terms of accuracy while increasing the width of the layer yields a significant increase in 
accuracy. 

As compared to the scheme put forward in [18] we observe a slight increase in the maximum error 
which is consistent with the observations made in [16] comparing the split and un-split PML methods for 
Maxwell’s equations. A direct comparison, however, is difficult due to differences in the computational 
scheme and boundary conditions. We emphasize, though, that the present results are arrived at without the 
use of filtering, thus confirming the stability of the scheme given in Eq.(13) and the associated analysis of 
well-posedness and decaying properties of the fields inside the layers. 

4.2. The Convecting Case. While the development of PML methods for the quiescent equations of 
acoustics relies on the analogy with the equations of electromagnetics, no such connection is possible in the 


13 








Fig. 2. u-velocity contours for a sound pulse , propagating in a quiescent medium . The contour levels are ±0.1, ±0.05, ±0.01 
and ±0.005 with the computed result given at t = 20 (a), t = 40 (b), t = 60 (c) and t = 80 (d). 


more general case of a convecting mean flow. 

The first idea that comes to mind is to introduce a new reference frame, moving with a speed of M 
along x and then apply the PML scheme developed in Sec. 4.1. This approach, however, has the unfortunate 
consequence that the layer interface becomes a moving interface in physical space. 

In [23] the use of a Lorentz-like transformation, connecting the convecting and quiescent wave-equations, 
is suggested in order to transform the quiescent PML method, such as to be useful in the convecting case. 
While this approach turns out to work well for the sound-waves, the resulting PML method has an abruptly 
changing convective velocity for the entropy and vorticity waves, resulting in significant reflection from such 
waves becoming quiescent exactly at the layer- interface. Moreover, the correct use of this approach in the 
corner regions of the PML layers is much less clear. 

Here we shall take a different approach although we shall rely on the PML schemes developed in the 
previous section combined with a few other techniques. Introducing layers in which the flow is accelerated into 
a supersonic region, thereby eliminating the need for absorbing boundary conditions, was recently proposed 
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Fig. 3. Maximum error at x = —48 as a function of time as computed with different types of boundary conditions and 
varying width of the PML layer. 

in [7] and modified in [9]. While this approach is appealing, it has an undesirable effect on the time-step of 
the whole computation and primitive sponge layers are still needed to yield an acceptable performance [9]. 

We propose to decelerate the flow, rather than accelerating it, to approach a quiescent flow inside the 
layer and then combine this approach with the PML scheme developed in Sec. 4.1. While such a scheme 
cannot be expected to be perfectly absorbing in the case of a finite layer, it does have the potential of 
a very efficient absorption, provided the deceleration is sufficiently smooth and the layer width is chosen 
accordingly. 

We propose to consider a PML- like scheme for the convecting case, Eqs.(l)-(2), of the form 


(23) 


Q t +M\l m(x)\ d P x - ^ ^ e'Q x p!Q y oMp , 


du r-, /vi 9u 

- + M\ l-m(x)]- 


dp 

- ~ 2 eu eP e , 

ox 


dv dv dp 

— +M[ 1 — m{x)\ — = - — - 2pv - p.P y - oMv , 
dp dp du dv 


■ ^ Qx M Qy i 


dP x 

dt 

dP, 


eu 


= pv 


dQx 

dt 

dQ y 


— &Qx w , 

= -flQy + V . 


dt ^ ’ dt 

Here e and // remains unchanged from Sec. 4.1. We have introduced m(x), which provides the decelerating 
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term by being m(a) — 0, with x — a signifies the layer interface, while m(a + x pm l) — 1 at the termination 
of the absorbing layer. We have found that using the error function provides a good compromise between 
steepness and smoothness such that 



where z = (x — o)/xpml and and x m controls the steepness and relative position of the profile, re- 
spectively. In Eq.(23) we have also introduced simple absorbing terms in the equations for p and v. Since 
the quiescent PML scheme only provides perfect absorption for the sound waves, this is meant to provide 
a simple mechanism for damping of the entropy and vorticity waves inside the layer. The parameters a(x) 
can be used to control the strength of this sponge layer for p and v. 

A few comments regarding the scheme, Eq.(23), is appropriate. First of all we note that for M = 0 we 
recover Eq.(13). Also since only the diagonal entries of A in Eqs. (l)-(2) are altered the well-posedness of the 
equations of acoustics remains intact. The philosophy here is that as the convective waves are slowed down, 
they approach the case of the quiescent acoustics for which Eq.(13) was shown to perform well. Moreover, 
slowly decelerating the waves as they enter the layer has the beneficial consequence that the wave fronts 
become increasingly normal to the boundary of the layer - much like water waves always being aligned with 
the beachfront. Hence, applying characteristic boundary conditions for truncating the PML layer can be 
expected to be efficient. 

4.2.1. A Numerical Example. In order to establish the soundness of the arguments that lead to the 
PML- like scheme given in Eq.(23), we have conducted a number of experiments using the scheme and initial 
conditions described in Sec. 4.1.1 and taking M = 0.5 as the convective Mach number of the mean flow. 

The decelerating term, Eq.(24) is generally specified by using <j m = 5 and x m — 0.5, i.e. the profile is 
centered in the middle of the absorbing layer. We have taken a(x) = e(:r), although this is by no means a 
unique choice and alternatives might well yield better performance that reported here. 

Since the layer now has multiple functions, i.e. it decelerates the waves while also acting as an absorbing 
layer, it is expected that, compared to the quiescent case, slightly wider layers should be used to achieve an 
acceptable accuracy. 

In Fig. 4 we show the temporal development of the density for the initial conditions given in Eq.(21) 
with e and p being given in Eq.(22) and the parameters chosen as in Sec. 4.1.1. We have taken the width 
of the layer as xpml = Vpml — 20, i.e. 20 computational cells, and At — 0.5. The exact solution is given 
in [22]. 

As expected, the sound pulse, as well as the entropy pulse, leaves the computational domain with no 
noticeable reflections from the layer. In Fig. 5 we show the development of the u-velocity component, 
arriving at similar conclusions. 

To address, in a more quantitative manner, the accuracy of the proposed scheme as a function of the 
width of the layer we have computed the maximum error in the pressure along the line x = 48 as a function 
of time. In Fig. 6 we plot the results for increasing width of the layer and compare to the accuracy obtained 
when using only a characteristic boundary condition to terminated the computational domain. 

Indeed we find that using a layer of only 10 cells yields an overall accuracy of the order of the approxi- 
mation error of the scheme and is superior to that obtained using characteristic boundary conditions only. 
By increasing the width of the layer to 20 cells, we observe a significant reduction, much like the case of the 
true PML in Fig. 3, of the reflections from the layer. 
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Fig. 4. Density contours for a sound and entropy pulse , propagating in a convecting medium with M = 0.5. The contour 
levels are ±0.1, ±0.05, ±0.01 and ±0.005 with the computed result given at t = 15 (a), t = 30 (b), t — 45 (c) and t = 60 (d). 


As expected, a slightly wider layer, as compared to the results in Sec. 4.1.1, is required in order to obtain 
an acceptable accuracy. However, rather than increasing the number of cells in the layer one could use a 
mapping, thereby stretching the grid, combined with a filter inside the layer. This approach was proposed 
in [6] for the case of acoustics and successfully used for the case of electromagnetics in [15, 16]. While this 
approach, certainly, will improve on the performance of the scheme with only little extra computational 
effort, we have chosen not to implement this technique in order to emphasize that the present schemes do 
not require the use of a filter in order to maintain stability. 

5* Concluding Remarks. The purpose of this paper has been two- fold. In the first part of the paper 
we provide an analysis of a recently proposed PML method for the equations of acoustics [18]. As remarked 
in [18, 19] these PML methods suffer from intrinsic numerical instabilities and we provided an explanation for 
this in terms of loss of well-posedness of the split set of equations and, as a result of this, the appearance of 
ill-posedness under small arbitrary perturbations. Such perturbations will inevitably exist in any numerical 
implementation of the split set of equations, rendering the schemes inherently unstable unless some kind of 
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Fig. 5. u-velocity contours for a sound and vorticity pulse, propagating in a convecting medium with. M = 0.5. The 
contour levels are ±0.1, ±0.05, ±0.01 and ±0.005 with the computed result given at t = 15 (a), t = 30 (b), t = 45 (c) and 
t = 60 (d). 


high-frequency damping, e.g. in the form of a low-pass filter, is introduced. 

The use of filters is a subject of some controversy. We believe, however, that while there might be 
numerous physical arguments for applying filters in various situations, it is somewhat of a concern if the 
numerical scheme, rather than the physics, dictates the need for a filter as is the case of the PML methods 
in [18]. Indeed, in situations where smooth initial conditions and only smooth boundaries are considered it 
is troublesome if the solution of a linear hyperbolic system requires the use of filters. 

In the second part of the paper we present a PML scheme for the quiescent equations of acoustics 
and prove that it is indeed absorbing for all frequencies and angles of incidence while remaining strongly 
well- posed. The properties of the layer is studied in more detail through numerical tests, confirming the 
analysis. 

In the general case of a convecting mean flow, the situation is no perfectly resolved. While we present an 
absorbing layer scheme that exhibit PML-like behavior, it is not a PML method, but rather a scheme arrived 
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Fig. 6. Maximum error at x = 48 as a function of time as computed with different types of boundary conditions and 
varying depth of the absorbing layer. 

at by combining several techniques. The scheme remains well- posed and performs well, although it requires 
the use of a slightly wider layer as compared to the true PML method presented earlier. The advantage of 
this scheme is that it contains the true PML scheme in the limit of a vanishing mean velocity and extends 
trivially to the general case of a mean flow which is not aligned with the axis. 

The development of a true well-posed PML method for the convective equations of acoustics remains a 
open challenge and we hope to address this in a future paper. 
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